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ABSTRACT 

We present a novel method for particle splitting in smoothed particle hydrodynam¬ 
ics simulations. Our method utilizes the Voronoi diagram for a given particle set to 
determine the position of fine daughter particles. We perform several test simulations 
to compare our method with a conventional splitting method in which the daughter 
particles are placed isotropically over the local smoothing length. We show that, with 
our method, the density deviation after splitting is reduced by a factor of about two 
compared with the conventional method. Splitting would smooth out the anisotropic 
density structure if the daughters are distributed isotropically, but our scheme allows 
the daughter particles to trace the original density distribution with length scales of 
the mean separation of their parent. We apply the particle splitting to simulations 
of the primordial gas cloud collapse. The thermal evolution is accurately followed to 
the hydrogen number density of 10^^ cm“^. With the effective mass resolution of 
^ 10“^ Mq after the multi-step particle splitting, the protostellar disk structure is 
well resolved. We conclude that the method offers an efficient way to simulate the 
evolution of an interstellar gas and the formation of stars. 

Key words: gravitation — hydrodynamics — methods: numerical — stars: formation 
— stars: Population III 


1 INTRODUCTION 

Smoothed particle hydrodynamics (SPH) is a widely used 
technique to follow the dynamics of a gravitationally inter¬ 
acting gas. A notable advantage of SPH is that the resolu¬ 
tion automatically increases with SPH particles concentrat¬ 
ing in dense regions. Nevertheless, in simulations of some 
hydrodynamical problems such as prestellar gas cloud col¬ 
lapse, one typically needs to follow a density evolution over 
too many orders of magnitude. In the central densest re¬ 
gion, SPH particles would e ventu ally violate a resolution re¬ 
quirement. iTruelove et al.l (|l997l) argue that a local Jeans 
length needs to be resolved by several grid cells at any 
time in a grid based hydrodynamics simulation. For a SPH 
simulation with the number of the neighbor particles iVngb 
and the particle mass rup, the resolution criterion may be 
written as Mmin < Mjeans/T^cr, where Mmin = AngbUlp 
is the minimal mass resolved by the SPH particles, and 
Mjeans ~ IS the Jcaus mass at the region 

with the sound speed Cs and the density p. Given a criti¬ 
cal mass ratio T^cr, one can determine the maximal particle 
mass required to resolve each simulated region. Insufficient 
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mass resolution results in producing artificial features such 
as smoothe d round shape, and of t en tri ggers unphysical frag¬ 
mentation (|Truelove et al.|[T997l . Il99^ i. Besides, it is costly 
to perform simulations with a uniform, and sufficiently high 
mass resolution so that the entire simulation region is fully 
resolved throughout the run. 


To save the computational time and memory, multi- 
step zoom-in techniques have been developed. Mesh re¬ 
finement can be done in a conceptually simple manner in 
grid-based simulations, where the cell division procedure 
can be systematically defined. In adaptive mesh refinement 
(AMR) simulations, for example, an initially large cell can 
be progressively refined with just smaller size rectangular 
cells. On the other hand, refinement of SPH particles is not 
trivial. A promising refinement technique for SPH simula¬ 
tions is particle splitting. A coarse parent particle which 
is about to violate certain resolution criteria is replaced 
with a number of finer daughter particles. A simple way of 
particle splitting is to distribute daughter particles within 
the smoothing length /iparent of the parent particl e (e.g. 
iKitsionas WhitworthI 20021: Bromm Loeij l2003l i. The 


splitting method of iKitsionas WhitworthI (|2002l i (h e 
after, KW02) is widely used i n the simulations of the super- 
massive blackhole formation (jPotti et al.ll200'^ : iMaver et al.l 





















2 G. Chiaki & N. Yoshida 


2 OIOII and of primo rdial star formation ([Yoshida et al.l[2006l : 
Hirano et al.|[20l3 h The KW02 method distributes thirteen 
daughters on a hexagonal closest packing array at a dis¬ 
tance /daughter = 1.5 /^parentfrom the Central particle. 
One of the daughters is located the same coordinate as the 
parent, and the angular orientation of the twelve daughters 
is randomly determined. Since there are, by definition, Nngh 
particles within /iparent, it can happen that a daughter parti¬ 
cle is placed very close to one of the neighboring particles by 
chance. This can result in significant overestimation of the 
local density. Also the fine anisotropic structures with sizes 
less than ~ /iparent tend to be smoothed out by imposing 
the spherical distribution on the daughter particles. Figure 
[ 1 ] illustrates an example particle distribution. The left panel 
shows the two-dimensional coordinates of the daughters gen¬ 
erated by the central coarse particle. After splitting with a 
random angle orientation, several daughters (red dots) hap¬ 
pen to be close to neighbors (black dots) within the smooth¬ 
ing kernel (orange-shaded region). In addition, daughters 
distributed in the originally underdense region can spuri¬ 
ously boost the local density estimate. 

The density perturbations would be mitigated or would 
cause only local effects if the daughter particles are placed 
in the region domina ted essentially by a single parent par¬ 
ticle. Along this idea, iMartel et al'l (|2Q06l i (MES06) present 
a method in which the daughter particles are placed on the 
vertices of the cube centered at a target parent particle with 
the cube side-length being the half inter-particle distance. 
We here explore yet another, and an elegant space parti¬ 
tioning method based on the Voronoi tessellation. A Voronoi 
cell is defined as a set of the points closest to each particle. 
The right panel of Figure [T] illustrates our scheme of parti¬ 
cle splitting in two dimension. The daughter particles (red 
dots) are distributed within a Voronoi cell (defined by the 
black lines). They are expected to reproduce the original 
density structure of the coarse parent particles. A nice fea¬ 
ture of the method is that the refinement can be done in 
a systematic manner for the given distribution of the orig¬ 
inal coarse particles. A number of advantages of using the 
Voronoi tessellation have been shown in p revious studies . 
For example, a moving mes h code AREPO ([SpringellbOlOl i 
and a mesh-free code GIZMO (|Hopkinsll2Q14l i implement the 
particle refinement based on the Voronoi tessellation. 

In the present paper, we first apply the particle split¬ 
ting based on the Voronoi diagrams to SPH simulations. We 
examine how our method improves the effective resolution 
while preserving original density structures, compared with 
a conventional splitting method of KW02 by several tests. 
As a realistic application, we perform simulations of gravi¬ 
tational collapse of a primordial gas cloud with two splitting 
methods. 


2 SPLITTING METHOD 

We construct the Voronoi diagram as follows. In our scheme, 
the Voronoi cell for each particle is obtained. We set a cube 
of size 2/cube centered at a target parent particle as a first 
trial Voronoi cell. The cube is then cut by a perpendicular 
bisector of a line segment which connects the central particle 
and each particle in the cube. After the cube is cut by all 
the bisectors, we obtain the Voronoi cell. Let us call this 
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Figure 1. We compare the isotropic particle splitting scheme 
(left) and our new scheme based on the Voronoi diagram (right) in 
two dimension. In the left panel, the central particle (blue point) 
is split to seven daughter particles (red points) distributed within 
the smoothing length of the parent particle (orange-shaded re¬ 
gion). In the right panel, the daughter particles are placed within 
the Voronoi cell of the central particle (blue-shaded region). 



Figure 2. Two-dimensional illustration of the coordinates of 
daughter particles (red dots) around the parent particle (blue 
dot). A daughter is on the center of mass of a subcell defined 
by a blue-dotted line which connects the parent and the middle 
point of each Voronoi edge. The plot is a zoom-in to the central 
portion of the right panel in Figure 1. 


procedure the facet-cut algorithm. Initially the cube size is 
set as /cube = /iparent- Particlcs outside the initial cube but 
closer to the central particle than 2/max can define the final 
shape of the Voronoi cell, where /max is the maximal distance 
between the central particle and the vertices of the Voronoi 
cell. If /max is larger than /cube/2, we enlarge the initial cube 
and repeat the facet-cut algorithm. This algorithm scales as 
0{Nsp log Nsp)-0{Nsp) for the number of split particles Agp, 
with a weak dependence on the distribution of the particles. 

There are various possibilities to determine the coor¬ 
dinates and mass of the daughter particles. Basically, the 
daughter particles of a single parent are desired to be well 
separated to each other and to have a nearly equal mass. 
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Figure 3. Top panels (a) and (b): Configuration of a subcell (or¬ 
ange region and blue-dotted lines) and the position of a daughter 
particle (red dot) around a parent particle (blue dot). Bottom 
panel (c): Stereograph of the daughter particles (red dots) replac¬ 
ing the central parent particle (blue dot) by our splitting method. 
Black solid lines indicate the Voronoi edges and grey dots indicate 
the neighbor particles around the target parent particle. 

Figure [2] shows an example set of the daughter particles. We 
split the Voronoi cell of the central parent particle (blue dot) 
into small subcells defined by the blue dotted lines. A sub¬ 
cell corresponds to each Voronoi vertex. A daughter is then 
placed at the center of mass of each subcell. In this way, we 
can realize an approximately equal mass distribution of the 
daughter particles. 

We extend the above procedure to three-dimensional 
cases. Figure [3] illustrates our procedure. We begin with the 
configuration given in the panel (a). We first draw a bound¬ 
aries of subcells on each Voronoi plane by the blue dotted 
line between the two points P and where P is the cen¬ 
ter of mass of the Voronoi plane, and E is the middle point 
of the Voronoi edge. We then define the base planes of the 
subcell corresponding to a vertex vi as the orange-colored 
region. The whole shape of the subcell is defined as the poly¬ 
hedron surrounded by the blue dotted lines in the panel (b). 
The polyhedron consists of the base planes with the thick 
blue-dotted lines in the panel (b), which is the same as the 
orange region in the panel (a), and the target parent parti¬ 
cle (blue dot in the panel (b)). Next, we merge the subcells 
whose corresponding vertices are connected by the shortest 
edges. We perform this merging procedure in order not to 
have daughter particles significantly close to each other. A 
Voronoi cell in 3D typically has 20-30 vertices. If the num¬ 
ber of the Voronoi vertices is more than ten, we repeatedly 
merge the closest subcells until the total number is reduced 
to ten. At this point we obtain a nearly uniform distribution 
of the daughters. We choose the threshold number A^th = 10, 
because we do not want to have a huge mass difference be¬ 
tween the parent and the daughters. Refinement by a factor 
of ten is a reasonable choice. We have actually performed an 
additional simulation in which subcells are not merged and 


Table 1. Density deviations just after the particle splitting for 
random particle distributions 


run 

Pmax 

Pmin 

P 

(7p 

UNIF-HALF-SPHERE 

3.76 

0.49 

1.12 

0.248 

UNIF-HALF-VORO 

1.65 

0.63 

1.04 

0.141 

UNIF-ALL-SPHERE 

4.03 

0.72 

1.17 

0.223 

UNIF-ALL-VORO 

1.89 

0.71 

1.07 

0.149 


have confirmed that the results are not significantly affected 
by the threshold number. We finally put a daughter particle 
at the center of mass of each (merged) subcell as shown by 
the red dot in the panel (b). 

The bottom panel in Figure [3] is a stereogram showing 
the distribution of the ten daughters (red dots) generated 
around the parent (blue dot). We equally divide the mass 
of the parent particle into rudaughter = '^^parent/10. We give 
the daughters the same temperature as that of the parent Q 
The velocity of the parent particle is simply assigned to the 
daughters so that the linear momentum and kinetic energy 
are strictly conserved. The daughters are given half the in¬ 
tegration timestep of their parent as a temporal value after 
splitting, but the individual integration timestep is recalcu¬ 
lated for the new set of particles in the next timestep. 


3 TEST SIMULATIONS 

We perform three test simulations to examine the splitting 
methods of the previous study (SPHERE) and ours (VORO). We 
are particularly interested in how the new density field after 
splitting differs from the original density field represented 
by the parent particl es. In the foll owing tests, we use the 
SPH code gadget-2 ([Springe II 2005 I) supplemented with the 
parallelized algorithms of the Voronoi tessellation and the 
particle splitting. All the control parameters are set to be 
the default values of gadget-2.0.7. To determine the SPH 
kernel size, we fix the number of the neighbor particles to 
E ngb = 50 zb 1. The densities of SPH particles are calculated 
in the standard manner using the cubic spline kernel. 

3.1 Random and Uniform Particle Distributions 

We first run simple simulations of an isothermal gas with¬ 
out gravity. We distribute 16^ SPH particles randomly in 
the cubic box of x G [0,1), y G [0,1), and z G [0,1) and 
adopt the periodic boundary conditions. The total mass 
of the particles is 1 and thus the mean volume-averaged 
density is pini = 1 in the code unit. The initial parti¬ 
cle mass is mp,ini = 1/16^. We set the isothermal sound 
speed ct such that the sound crossing time is equal to 
the unit time. Since the typical length scale of the den¬ 
sity fluctuation resulting from the random noise is compa¬ 
rable to the smoothing length of the particles, we define the 
characteristic sound crossing time as tsc = ^ini/cT, where 

^ We have a few other choices such as assigning entropy. We 
choose temperature as a primary thermodynamic quantity for 
simplicity. 
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hint = {3Nnsh'mp,ini/4:TrpinG^^ is the initial mean smoothing 
length. We first run a simulation to t = 10 without particle 
splitting to reduce the random noise by the gas pressure. At 
the time, the deviation of the density <jp is reduced down to 
0.0026. The maximal and minimal densities are found to be 
(pmax,Pmin) = (1.08, 0.92). Wc Call thc sct of the runs UNIF. 

We then split half of the particles in x G [0,0.5) for 
run HALF and the all particles for run ALL. Each set is per¬ 
formed with the two splitting methods of SPHERE and VORO. 
After particle splitting, the root-mean-square of the density 
(Tp suddenly increases; the maximal density increases and 
also the minimal density decreases. Note that the particle- 
averaged density fluctuates around unity after perturbed by 
the particle splitting. The density deviation eventually ap¬ 
proaches to zero and the minimal and maximal densities 
approach to the mean density because of the gas pressure. 
We summarize the results in Table [H The second to fifth 
columns show the maximal, minimal, average, and root- 
mean-square of the density just after the particle splitting, 
respectively. 

For the run HALF, the density deviation from the initial 
value is smaller with VORO than with SPHERE. The maximal 
and minimal densities just after the splitting change from 
(pmax,Pmin) = (3.76,0.49) to (1.65,0.63), with the latter 
being closer to 1. The density deviation is also reduced from 
0.25 to 0.14. The average density is larger than the original 
density by 12 % for SPHERE, while 4 % for VORO. In the region 
X G (0.1, 0.4) populated with the fine particles, the density is 
overestimated at most by a factor of 4 and underestimated 
by a factor of 0.3 for SPHERE, whereas the density is overes¬ 
timated by a factor of 0.6 and underestimated by a factor 
of 0.4 for VORO. At the boundary of the split and un-split 
particles, the density of the coarse particles is overestimated 
by a factor of 2.5 and underestimated by a factor of 0.5 for 
SPHERE, while the density deviates by a factor of 0.3 for VORO. 
At the time t = 20, when we stop the simulations, the den¬ 
sity difference still remains at the boundary. The maximal 
and minimal densities are 0.3 and 0.2 for SPHERE and VORO, 
respectively. Because the particles with different masses are 
mixed near the boundary, the local density does not quickly 
converge to the mean value. In practical gas collapse simula¬ 
tions, however, the densest cloud core shrinks more rapidly 
than the boundaries of fine and coarse particles. Thus the re¬ 
maining density perturbation on the boundaries would have 
little effect on the central regions. 

For the run ALL, the original density field is well pre¬ 
served by our splitting method. The maximal density is 1.89 
for VORO, which is by a factor of two better than SPHERE. Al¬ 
though the density is slightly more underestimated for VORO. 
the overestimation of the average density is improved from 
17 % to only 7 %, and the root-mean-square of the density is 
reduced from 0.22 to 0.15. As we have seen in Section [U the 
density is overestimated for a particle that is placed signif¬ 
icantly close to one of the neighbor particles. This happens 
with the SPHERE configuration because it involves random¬ 
ization of the orientation of the daughter particle distribu¬ 
tion. With VORO, such an unfortunate case with significant 
density overestimation does not occur because the coordi¬ 
nates of the daughters are uniquely determined with the 
separation of any two daughters being a fraction of the dis¬ 
tance of their parents. We continue the run and find that 
another 4.6 sound-crossing time is required for the density 



Radius [pc] 

Figure 4. Density profile as a function of the radius for our test 
GRAD. The black dot-dashed curve shows the density profile before 
the particle splitting. The solid and dashed lines show the density 
profiles just after the particle splitting with the methods of the 
previous (red dashed) and this (blue solid) works, respectively. 

deviation to recover to the initial value 0.026 for both SPHERE 
and VORO. 

3.2 Particle Distribution with Density Gradient 

It is important to test our method in the case with a steep 
density gradient. If the length scale of a density fluctuation 
Ip = pl\dpldr\ is smaller than hparent ~ (A^ngb^p/p)^^^, the 
structure with length scales < /iparent would be smoothed 
out with the previous splitting method. On the other hand, 
with our method, the length scale over which the daughter 
particles are distributed is comparable to the inter-particle 
distance Ap ~ (mp/p)^/^ before split. Thus, the length over 
which the density profile is affected is expected to be much 
smaller. To illustrate this, let us consider a spherical cloud 
with a large density gradient p oc corresponding to 

small Ip. Such structure can be seen on the edge of the disk 
in rotating isothermal clouds fSection 13.3D . and also in a 
volume within a primordial gas cloud where the transition 
lines of HD molecules becomes optically thick (see Section 
|4|. We call this test GRAD. 

We realize the desired density profile with p oc by 
perturbing particles that are homogeneously distributed ini¬ 
tially. The black dot-dashed curve in Figure |4] shows the 
generated density profile. The sound speed of the cloud is 
Cs = 1 km s“^, corresponding to T = 100 K, actually be¬ 
ing close to the primordial gas cloud temperature. Then, 
we split the particles with densities nn > 3 x 10® cm“^. 
At the density, the length scale of the density variation is 
Ip = r/4 = 0.006 pc. The red dashed curves and blue solid 
curves in Figure |4] show the density profiles just after the 
particle splitting of the previous and our works. The den¬ 
sity profile is flattened at r = 0.03 pc over the length scale 
~ ^parent = 0.01 pc for SPHERE, whilc we find only a small 
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Figure 5. Density slice on the equatorial plane for the Boden- 
heimer test simulations without particle splitting (NoSp: top) and 
with the particle splitting of the previous work (SPHERE: bottom 
left) and the present work (VORO: bottom right) when the maximal 
density is ~ 5 x 10^^ cm“^. 


Table 2. Density deviations for BB79 tests 


run 

Np,ini 

Np^fin 

Pmax 

P 

Pmin 

P 

gp 

p 

BB79-NoSp 

BB79-SPHERE 

113,098,912 

4,188,896 

15,925,664 

5.54 

0.41 

0.33 

BB79-V0R0 

4,188,896 

12,805,187 

2.22 

0.62 

0.15 


bump with ~ Ap ~ 0.004 pc with our method. The daughter 
particles located within the inter-particle distance of their 
parent can preserve the density structure before splitting. 
Note that we consider a spherical cloud in this test. For 
a cloud with anisotropic density structures, distributing the 
daughters isotropically results in smoothing the original den¬ 
sity structure. The artificial smoothing is mitigated with our 
method, as we will see in the case of a primordial cloud in 
Section [H 

3.3 Bodenheimer test 

Next, we perform a we ll-known test introduced by 
IBoss Bodenheimerl (|l979l i (BB79). We follow gravitational 
collapse of an isothermal cloud with and without the parti¬ 
cle splitting. The cloud mass is Mini = 1 M© and its radius 
is Rini = 0.02 pc initially. The isothermal sound speed of the 
gas is ct = 0.17 km s“^. The ratio of the thermal energy to 
the gravitational energy is aini = 5clRini/2GMini = 0.32. 
We set the initial angular velocity flini = 7.2 x 10“'^® s“^ and 
then the ratio of the rotational energy to the gravitational 
energy is ^^ini = = 0.28. The amplitude of 

m = 2 perturbation is Am=2,ini = 0.1. The particles are dis- 



Figure 6. Radial profile of the density for Bodenheimer test sim¬ 
ulations without the particle splitting (NoSp: black dot-dashed) 
and with the particle splitting of the previous work (SPHERE: red 
dashed) and the present work (VORO: blue solid). Note that the 
red and blue curves closely lie on each other in the plot. 

tributed initially on a regular grid, and the particle mass is 
modulated as ?Tip,ini = 'mp,ini(l + ^m=2,ini cos 20/2), where 
^p,ini is the initial average particle mass, and 0 is the az¬ 
imuthal angle of the particles. 

Th e semi-analytic approaches by iTsuribe fc Inutsri^ 
(|l999bl i show that, in the cloud with (a,/3) = (0.32,0.28), 
the non-axisymmetric structure of the gas evolves, and the 
fragmentation pro perty depends on th e initial amplitude of 
the perturbation. iTsuribe Inutsii^ (|l999al i perform the 
three-dimensional simulations and conclude that such cloud 
fragments if the initial amplitude is larger than 0.15. In the 
BB79 test, we run simulations with the minimal mass reso¬ 
lutions IZcT = 1, 10, and 100. For IZcr = 1 and 10, artificial 
fragmenta tion occurs for the la ck of the resolution as sug¬ 
gested by iTruelove et al ] (Il997f) . We here show the results 
for IZcT = 100. 

We perform three runs and compare them. In the run 
NoSp, we distribute a very large number of particles (iVp ~ 
10®) to satisfy the Jeans criterion 71 = Mjeans/Mmin ^ 100 
until the gas density reaches 5 x 10“^^ g cm“®, where the 
gas would become optically thick. In the runs SPHERE and 
VORO, we initially distribute A^p,ini ~ 4 x 10® particles. When 
particles meet the splitting criteria, i.e., violating the Jeans 
condition, they are split one by one in the course of the 
simulations. 

Figure [5] shows the density slice of the clouds for NoSp 
(top), SPHERE (bottom left), and VORO (bottom right) at the 
maximal density pcen = 5 x 10“^^ g cm“®. The stretched 
thin filament does not fragment in al l the three runs up 
to the last output time as reported bv ITsuribe Inutsukal 
(|l999bl i. Figure [6] shows the evolution of the density profiles 
as a function of the distance from the densest part of one 
end of the filament. The overall structure of the filament 
remains essentially the same among the three runs. We con¬ 
clude that the resolution Tier = 100 is sufficient when we 
consider collapse and fragmentation of the cloud. The num¬ 
ber of particles is A^p,fin ~ 10^ for SPHERE and VORO when we 
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terminate the simulations at pcen = 5x 10“^^ g cm“®, which 
is by an order of magnitude less than NoSp. This means that 
the computational cost is well reduced by the splitting tech¬ 
nique although the final snapshots are similar to the one for 
NoSp. 

Let us compare the two runs SPHERE and VORO in de¬ 
tail. We follow the density deviations of the 1000 parti¬ 
cles that are split first. The results are summarized in Ta¬ 
ble [21 The particles are split for the first time when p = 
3 X 10~^^ g cm“^. Just before the splitting, the maximal and 
minimal densities relative to the mean density of the 1000 
particles are 1.08 and 0.85, respectively, and the deviation of 
the density is 0.037. Just after the 1000 particles are split, 
the density deviation marks the maximum. In the SPHERE 
and VORO runs, the maximal and minimal densities just af¬ 
ter the particle splitting are (pmax/p, pmin/p) = (5.54,0.41) 
and (2.22, 0.62), respectively. The maximal deviations of the 
density from the original value 1 are reduced more than by 
a factor of two by our splitting method. The density devia¬ 
tion relative to the mean density just after the splitting also 
decreases from 0.33 (SPHERE) to 0.15 (VORO), approaching 
to the initial value of 0.037. The result shows that, also in 
the practical collapse simulation, the daughter particles dis¬ 
tributed in a Voronoi cell of each split particle are not close 
to other neighbor particles, and their density is not signif¬ 
icantly overestimated, while this happens with the conven¬ 
tional splitting method. 

In this BB79 test, the density structure is not very differ¬ 
ent between SPHERE and VORO, showing an interesting con¬ 
trast to our test GRAD that is more relevant to a primordial 
warm gas (see Section 13.21) . In the cold gas with small ct 
in BB79, the smoothing length hparent ~ 
is smaller than the length scale Ip all over the region if the 
critical condition is set as IZcv = 100. On the iso-density 
surface of p = 5 x 10“^^ g cm“^, where the second particle 
splitting occurs. Ip is smallest, 4 x 10“® pc, at the edge of 
the iso-density filament. The smoothing length in the region 
is /^parent = 3 X 10“® pc. SPH particlcs which satisfy the 
Jeans criterion IZcr = 100 can resolve the density structure 
when the gas is sufficiently cold (ct ~ 0.1 km s“^). This is, 
however, not in the case of the primordial warm gas as we 
revisit with detailed discussion later in Section 


4 COLLAPSE SIMULATIONS OF 
PRIMORDIAL GAS CLOUD 

Finally, we apply our splitting method to a more realistic 
calculation of the collapse of a primordial minihalo. We iso¬ 
late the spherical region with length ~ 1 kpc centered at 
one of the minihalos for med at the redshift z = 15 in the 
cosmological simulation of lHirano et al.l (|2014| j . The mass of 
the halo is Mh = 4.8 x 10^ M©. The halo is categorized to a 
subclass of the primordial minihalos which hosts less mas¬ 
sive stars (10-100 M©) because of the efficient cooling by HD 
molecules. The chemic al reaction ra t es an d radiative cool¬ 
ing rates are taken from I Chiaki et a DioH) for a primordial 
composition. We reduce the cooling rates of emission lines 
by t he Sobolev length method in the optically thick regime 
as in iHirano Yoshidal (|2013l b The splitting methods of the 
privious work (SPHERE) and this work (VORO) are employed 
for the gas particles. The number of neighbor particles is 



Figure 7. Distribution of the particles on the thermodynamic 
phase plane in the simulations of primordial gas cloud formaton, 
when the central density nH,cen = 10^ cm“^. The color indicates 
the mass-weighted number of the particles. Left and right panel 
show the result with the particle splitting of SPHERE and VORO, 
respectively. The dotted lines attached to the orders ffith’ indi¬ 
cate the density and temperature where the particles are split for 
the nth time. The SPH particles in the region below the white 
diagonal line have temperatures that are underestimated due to 
the particle spitting. 



Figure 8. Projected density distribution of the primordial clouds 
from the face-on and edge-on views of the disk methods of SPHERE 
(top) and VORO (bottom) when the central density is 10^^ cm“^. 


restricted to be Angb = 64 ±8, and the criterion for the par - 
ticle splitting is IZcr = 1000 as set bv iHirano et al ] il2ni4i . 
The central particles are already once split by SPHERE split¬ 
ting when the zoomed initial conditions are generated. The 
initial particle mass is = 5/13 M©. The number of 

gas particles is initially 169,511. When the central density 
reaches nn = 10^^ cm“^, the particles have been split at 
most five times and the number of split particles increases 
up to 314,867 for SPHERE and 318,148 for VORO. 

Figure [7] shows that distribution of the gas particles in 
the density-temperature plane. We use the snapshot when 
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Figure 9. Radial density profile of the primordial cloud with the 
particle splitting methods SPHERE (red dashed) and VORO (blue 
solid). We also plot the grey dotted line oc which a primor¬ 

dial gas cloud tends to be along. 


the central density nH,cen = 10^ cm“^. Note that HD 
cooling operates at the densities nu > 10^cm“^, where a 
disk structure is formed with noticeable spiral arms. The 
particles are split at the densities and temperatures indi¬ 
cated by the dotted lines in Figure [71 The index ‘nth’ at¬ 
tached to the dotted lines indicates the criterion for the 
nth particle splitting: Mjeans = 7^criVngb(^sPH/13’^~^) and 
'R'crNngh{m^p^/10^~^) for SPHERE and VORO, respectively. 
We see that the temperature significantly fluctuates along 
the adiabatic line T oc from the splitting points for 

SPHERE. Since the gas cooling and heating do not work 
promptly within a short timestep, the temperatures of 
daughters deviate along the adiabatic line when their den¬ 
sities are perturbed from the parent’s density one timestep 
after splitting. Thus the temperature deviation almost di¬ 
rectly reflects the density perturbation induced by the par¬ 
ticle splitting. The white lines in Figure [7| indicate the ex¬ 
pected trajectories, i.e., the gas particles having the lowest 
temperature at each density would evolve roughly on the 
white line without the particle splitting. The figure shows 
that our splitting method slightly improves the temperature 
deviation from the splitting point. The distributions of par¬ 
ticles on the density-temperature planes are not signihcantly 
different for two splitting methods SPHERE and VORO. 

We hnd that substantial differences are caused in the 
global structure by the two different splitting methods. Fig¬ 
ure [8] shows the snapshots of the gas cloud at the central 
density nH,cen = 10^^ cm“^. Comparing the middle panels, 
we hnd that the spiral arms have more developed in VORO 
than in SPHERE. This is related to the associated smoothing 
and isotropization after split, as we discussed in Section [321 
At densities nu ~ 10^-10^ cm“^, where the emission lines 
from HD molecules becomes optically thick. Ip is reduced to 
~ 0.01 pc because of the rapidly increasing pressure. The 
length at which the daughter particles are distributed is 
~ /^parent ~ (A^ngb^p/p)— 0.01 pc iu ruu SPHERE and 
~ Ap ~ (rup/py^^ ~ 0.004 pc in run VORO. Even for the 
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Figure 10. Distribution of the SPH particles above (orange) and 
below (blue) the white line in Figure [7| on the iso-density sphere 
with nil = 10^ cm“^ for the splitting methods of the previous 
(top) and this work (bottom). 


same Jeans criterion, our splitting method allows to make 
the resolution of daughter particles effectively higher than 
in SPHERE by a factor of ~ Aparent /Ap ~ ■ 

Figure [9] shows the radial density prohle and its evolu¬ 
tion. Our splitting method maintain the steep density prohle 
at the densities nu ~ 10® cm“^ where the emission lines of 
HD molecules are optically thick, and at ~ 10® cm“® where 
the gas heating by H 2 molecular formation via the three- 
body reactions is dominant. The right panels of Figure [8] 
show that the bar-like structure has developed in VORO while 
the central region appears smooth and spherical in SPHERE. 
Clearly, details of splitting method affect the evolution of 
the global structure such as spiral arms and also the frag¬ 
mentation of the cloud core. 

We suspect from the tests GRAD and BB79 that the 
anisotropic density structure in the primordial collapsing 
gas is smoothed by the isotropic distribution of daughter 
particles. In a warm gas cloud with large ct such as a pri¬ 
mordial gas cloud, the length Ip can not be resolved by 
^parent cx: (G^p) ~cvcu though the Jeans criterion 

with IZcT = 1000 is satished. To avoid the isotropic diffusion 
of density structures in a warm gas where ct ^ 1 km s“^, 
we should impose a severer Jeans criterion with larger res¬ 
olution threshold IZcr such that Ip can be resolved by the 
coarse particles with the conventional splitting method (see 
Section[5]). Our method overcomes this problem by distribut¬ 
ing the hne particles within an inter-particle distance from 
the parent, which shortens the effective resolution of the gas 
particles. Therefore, it is not necessary to have severer cri¬ 
teria for splitting than the Jeans criterion. 
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Figure 11. Distribution of the Toomre Q parameter in the disk 
from the face-on view for runs SPHERE (left), CUBE (middle) and 
VORO (right) in the same region as the middle panels of Figure [S] 
The Q parameter is smallest at the white squares in the spiral 
arms. 


5 DISCUSSION 

We have developed a novel method for particle splitting in 
SPH simulations based on the Voronoi diagram (VORO). Our 
series of tests show that the new method preserves well the 
original density structure before splitting, while achieving 
an effectively high resolution. The density perturbations in¬ 
duced by the splitting are reduced by a factor of two com¬ 
pared with the conventional splitting method (SPHERE). We 
have found that similar improvements are achieved by an¬ 
other particle splitting method of MES06 (CUBE hereafter), 
where daughters are distributed within a distance /iparent /8 
(^ Ap/4) from their parent. In the test GRAD, the density 
profile is similar to VORO. Only the slight bump appears at 
the splitting point over the length scale ~ Ap. In the test 
BB79, the morphology and the density profile is similar to 
runs SPHERE and VORO. The density deviation just after the 
splitting is better than VORO. The maximal and minimal den¬ 
sities relative to the average density of 1000 particles first 
split are 1.91 and 0.68, respectively. The root-mean-square 
of the density is 1.03. 

In the case of a contracting primordial gas cloud, where 
the gas no longer evolves isothermally, the temperature de¬ 
viates along with the density perturbation generated by the 
particle splitting as Figure [7] shows. Some particles are over¬ 
cooled and might trigger the artificial fragmentation by their 
small pressure, but this would not be in the case. Figure fTOl 
shows the spatial distribution of the SPH particles on the 
iso-density surface with nn = 10^ cm“^. We divide the gas 
particles into hot (orange) and cool (blue) components re¬ 
spectively above and under the white line in Figure [T] The 
blue dots indicate the overcooled gas. If the overcooled gas 
particles cluster, they might trigger the fragmentation by 
their small pressure. However, the cool gas particles appar¬ 
ently scatter in all directions for SPHERE (upper panel). For 
VORO, although the cool component appears to cluster at 
longitude 180°, where a spiral arm lies, the hot gas particles 
are also in the arm region and would prevent the collapse 
of the cool component. The local perturbations of temper¬ 
ature and density triggered by the particle splitting would 
not affect the gas fragmentation irrespective of the splitting 
methods. 

On the other hand, the global density structure such 
as spiral arms and bars can affect the gravitational stabil¬ 
ity of the gas. As seen in Figure (8] the high-density region 
remains strongly elongated in VORO run. Although the long¬ 
time evolution of the spiral arms is not followed in our sim¬ 


ulations, we can derive a reasonable estimate on whether 
or not the disk is stable by calculating the Toomre Q pa¬ 
rameter Q = 2ncT where S is the column density in 

each part on the disk. Figure ITT] shows the distribution of Q 
parameter. Clearly, the value of Q is small along the spiral 
arms and is the smallest at the point indicated by the white 
square in each panel: Qmin = 0.59 for SPHERE and 0.41 for 
VORO. We have also performed the collapse simulation with 
the splitting method CUBE. The distribution of Q parameter 
is shown in the middle panel of Figure [11] The elongation 
of the central core is slightly smaller than our method VORO. 
The inner structure of the spiral arms is resolved as well as 
VORO, but the density structure appears slightly more dif¬ 
fuse. This is likely because of the nature of the method of 
CUBE, where the daughter particles are distributed on the 
vertices of the cube around a parent, and the original density 
structure formed by parent particles are slightly smoothed 
in the directions both perpendicular and parallel to the disk. 
The column density in the spiral arms for the run CUBE is 
less than VORO, and the resulting minimal Q parameter is 
0.58. The l inear analysis of the ring-mode instabilities by 
iTakahashU (|2015l i suggests that spiral arms with Q < I/tv 
is dynamically unstable and likely produce fragments. We 
can say that the spiral arms are marginally stable against 
fragmentation in the both runs at the time of the snapshot. 
Still, the difference in the global structure of the protostellar 
disk affects the evolution of the central protostar in a com¬ 
plica ted manner (see e.g. iGreif et 'n]l2012l : I Vorobyov et al.l 
120131 ) . It is important to fully resolve the disk structure over 
a long time. 

Our test simulations suggest that the mass resolution 
with VORO is effectively higher than with SPHERE by a factor 
of ~ /iparent /Ap ~ ^ngb' small-scalc auisotropic 

structures are smoothed out when the particles are split in 
the region with steep density gradient with insufficient res¬ 
olution with SPHERE method. To avoid the artificial mix¬ 
ing, we suggest to impose an additional criterion than the 
Jeans criterion for splitting with the conventional splitting 
method with spherical distribution of daughters: the length 
scale of the density fluctuation Ip should be sufficiently re¬ 
solved by the local SPH smoothing length. To confirm this, 
we have performed another simulation with SPHERE splitting 
with a substantially higher resolution of IZcv = 4000. The 
high-resolution run indeed shows an elongated cloud and 
apparent spiral arms, similarly to those found in the VORO 
run. We emphasize that our splitting method performs only 
with the modest Jeans criterion, because the distribution 
of daughters traces the original anisotropic density struc¬ 
ture before splitting. The method saves computational cost 
by effectively reducing the number of particles necessary to 
resolve the cloud structures. 

It is actually costly to construct the Voronoi diagrams 
with the current algorithm. In the simulations with IZcr — 
1000, the computational time until the density of the cloud 
center reaches nH,cen = 1 x 10^^ cm“^ is 48,315 (SPHERE) and 
81,320 (VORO) seconds using 128 cores. The Voronoi tessel¬ 
lation accounts for about a half of the total computational 
time. However, in the simulation of SPHERE with IZcr = 4000, 
where the bar-like structure is seen also with the conven¬ 
tional splitting method, the computational time is 323,309 
seconds, exceeding the time for VORO with IZcr = 1000. 
Clearly, the formation of small-scale structrue such as the 
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bar-like structure is accurately followed with smaller com¬ 
putational time with our splitting method. We can further 
save the computational cost for generation of the Voronoi 
tessellation by introducing the i ncremental construct ion al¬ 
gorithm scaling as 0{N log N) (|Sugihara Iril[l994l l . In a 
forthcoming paper, we present an efficient scheme to paral¬ 
lelize the algorithm that scales as 0{N) (Chiaki, Yoshikawa 
V Yoshida in prep.). 

We note that there have been some simplifications 
adopted in our method. First, we place the daughter par¬ 
ticles at the center of the mass of the subcells as shown in 
Figure [3l Another promising choice might be to place daugh¬ 
ters such that the original Voronoi planes between the parent 
and the neighboring particles are not affected even after the 
splitting. Second, we still follow th e traditional SPH s cheme 
to solve the equation of motion. ISaitoh fc Making (|2013l i 
point out that the traditional SPH scheme causes asymmet¬ 
ric errors of density and pressure estimations at boundaries 
between the two phases with different densities. The density 
and pressure errors obviously affect the fluid motion. We find 
that, in the case of our particle splitting, the pressure pertur¬ 
bation at the boundary between the fine and coarse regions 
nearly cancel out and thus spurious force is not induced. 
Third, the distribution of the parent particles has intrinsic 
noise due to its discrete nature. Then unphysical anisotropic 
shapes of the Voronoi cells might generate and even pro¬ 
mote density perturbations. We could re gulate the sh apes of 
the Voronoi cells by Lloyd’s algorithm (|Llovdlll982l i. which 
is expected to yield accurate density estimates. Finally, we 
adopt the cubic spline kernel to estimate the density as in 
the standard SPH. The choice of other kernel functions and 
the number of the neighbor particles would affect the accu¬ 
racy of the density reconstruction. Overall, we have shown 
that our splitting method based on the Voronoi diagram im¬ 
proves the preservation of the density and the morphology 
of gravitationally contracting clouds even with the minimal 
implementation, but further studies are certainly warranted 
to optimize the choice and the implementation of each part 
of the splitting scheme. 

The splitting method can be applied not only to simu¬ 
lations of star forming clouds but also to other astrophysi- 
cal problems. For example, in the downstream of the cool¬ 
ing shocks, th e density rapidly increases by three orders 
of magnitude (jChiaki et al.ll2013 L The density structure is 
highly anisotropic around the shock front by hydrodynamic 
instabilities. The dense shocked region can be fully resolved 
by our splitting method while maintaining the fine struc¬ 
ture. We can also perform de-refinement of the particles as 
the converse operation of the particle splitting. The parti¬ 
cles to be merged can be uniquely identified once we de¬ 
fine the Voronoi diagram. An exampl e implemented in a 
moving-mesh code is found in AREPO (|Springelll2010l L The 
de-refinement can also be applied to the problem of an inter¬ 
stellar shock, where the density of fluid elements decreases 
after passing through the shocked region. In such cases, de¬ 
refinement of particles behind the shock will significantly 
save the number of particles and the overall computational 
time. We continue developing the de-refinement techniques 
based on the Voronoi diagrams. Our method of particle spit¬ 
ting presented here and de-refinement technique can be ap¬ 
plied to particle-based simulations of various problems as 


elaborated schemes both to preserve complex density struc¬ 
tures and to save the computational cost. 
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